library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
library(sdmTMB)
library(tidyr)
Mesh
# 300 knots
plot(mesh)
Model specifications
m1 <- sdmTMB(
bio.adult ~ 0 + as.factor(Year) + GearCat + s(Month, bs = "cc", k = 10),
data = dat,
mesh = mesh,
offset = log(dat$SweptArea),
family = tweedie(link = "log"),
spatial = "on",
time = "Year",
spatiotemporal = "IID"
)
## Warning in checkDepPackageVersion(dep_pkg = "TMB"): Package version inconsistency detected.
## glmmTMB was built with TMB version 1.8.1
## Current TMB version is 1.9.0
## Please re-install glmmTMB from source or restore original 'TMB' package (see '?reinstalling' for more information)
## Warning: The model may not have converged. Maximum final gradient:
## 0.0108123051405755.
Some diagnostics
dat$resids <- residuals(m1) # randomized quantile residuals
hist(dat$resids)
qqnorm(dat$resids[is.finite(dat$resids)])
qqline(dat$resids[is.finite(dat$resids)])
# check spatial autocorrelation
ggplot(filter(dat,Year%%5==1), aes(lon, lat, col = resids)) +
scale_colour_gradient2() +
geom_point() + geom_sf(data=map,inherit.aes = FALSE)+xlim(-12,30)+ylim(47,62)+
facet_wrap(~Year,ncol=2)
## Warning: Removed 758 rows containing missing values (geom_point).
Terms
m1
## Spatiotemporal model fit by ML ['sdmTMB']
## Formula: bio.adult ~ 0 + as.factor(Year) + GearCat + s(Month, bs = "cc",
## Formula: k = 10)
## Mesh: mesh
## Time column: Year
## Data: dat
## Family: tweedie(link = 'log')
##
## coef.est coef.se
## as.factor(Year)1970 11.34 1.24
## as.factor(Year)1971 10.15 0.69
## as.factor(Year)1972 10.10 0.67
## as.factor(Year)1973 10.01 0.83
## as.factor(Year)1974 9.99 1.43
## as.factor(Year)1975 9.08 0.82
## as.factor(Year)1976 10.04 0.74
## as.factor(Year)1977 9.53 0.71
## as.factor(Year)1978 11.58 0.52
## as.factor(Year)1979 11.61 0.51
## as.factor(Year)1980 12.24 0.45
## as.factor(Year)1981 13.35 0.47
## as.factor(Year)1982 11.87 0.46
## as.factor(Year)1983 12.14 0.45
## as.factor(Year)1984 12.45 0.44
## as.factor(Year)1985 11.03 0.43
## as.factor(Year)1986 11.23 0.42
## as.factor(Year)1987 10.88 0.42
## as.factor(Year)1988 11.42 0.41
## as.factor(Year)1989 12.21 0.41
## as.factor(Year)1990 10.85 0.40
## as.factor(Year)1991 10.92 0.40
## as.factor(Year)1992 10.67 0.40
## as.factor(Year)1993 11.15 0.39
## as.factor(Year)1994 10.88 0.38
## as.factor(Year)1995 10.50 0.38
## as.factor(Year)1996 10.38 0.38
## as.factor(Year)1997 10.63 0.37
## as.factor(Year)1998 10.63 0.38
## as.factor(Year)1999 10.60 0.37
## as.factor(Year)2000 10.32 0.37
## as.factor(Year)2001 10.08 0.36
## as.factor(Year)2002 9.61 0.37
## as.factor(Year)2003 9.65 0.35
## as.factor(Year)2004 9.71 0.36
## as.factor(Year)2005 9.46 0.36
## as.factor(Year)2006 9.57 0.36
## as.factor(Year)2007 9.93 0.36
## as.factor(Year)2008 10.16 0.35
## as.factor(Year)2009 9.47 0.36
## as.factor(Year)2010 9.92 0.35
## as.factor(Year)2011 10.17 0.34
## as.factor(Year)2012 9.84 0.35
## as.factor(Year)2013 9.72 0.35
## as.factor(Year)2014 9.80 0.35
## as.factor(Year)2015 9.81 0.35
## as.factor(Year)2016 9.93 0.35
## as.factor(Year)2017 9.59 0.35
## as.factor(Year)2018 9.16 0.36
## as.factor(Year)2019 8.59 0.36
## as.factor(Year)2020 8.47 0.37
## as.factor(Year)2021 8.32 0.36
## GearCatBT 0.25 0.18
## GearCatGOV_CL 1.54 0.17
## GearCatGOV_GG 0.97 0.22
## GearCatOTT -3.26 0.53
## GearCatPEL -2.05 0.75
## GearCatTV 0.89 0.23
##
## Smooth terms:
## Std. Dev.
## sds(Month) 0.43
##
## Dispersion parameter: 665.08
## Tweedie p: 1.61
## Matern range: 1.35
## Spatial SD: 0.37
## Spatiotemporal SD: 0.26
## ML criterion at convergence: 248325.051
##
## See ?tidy.sdmTMB to extract these values as a data frame.
Predict onto the grid
grid=grid %>% filter(!Year %in% c(1967,1968))
grid$Year=as.integer(grid$Year)
dat %>% group_by(GearCat) %>% summarise(n=n())
## # A tibble: 7 × 2
## GearCat n
## <fct> <int>
## 1 BAK_GG 520
## 2 BT 33421
## 3 GOV_CL 37009
## 4 GOV_GG 6353
## 5 OTT 87
## 6 PEL 88
## 7 TV 10896
grid$GearCat = 'GOV_CL'
# calculate the mode of month
# getmode=function(x){uniq=unique(x)
# uniq[which.max(tabulate(match(x,uniq)))]}
# getmode(dat$Month) Month set to 9
grid$Month = 9
p <- predict(m1, newdata = grid)
Spatial random fields
plot_map(p, "omega_s") +
scale_fill_gradient2() +
ggtitle("Spatial random effects only")
## Warning: Removed 4817 rows containing missing values (geom_raster).
Spatio-temporal random fields
plot_map(filter(p,Year%%5==1), "epsilon_st") +
scale_fill_gradient2() +
facet_wrap(~Year,ncol = 2) +
ggtitle("Spatiotemporal random effects only")
## Warning: Removed 1169 rows containing missing values (geom_raster).
plot_map(filter(p,Year%%5==1), "exp(est)/1000") +
scale_fill_viridis_c(
trans = "pseudo_log",
# trim extreme high values to make spatial variation more visible
na.value = "yellow", limits = c(0, quantile(exp(p$est)/1000, 0.995)),'cod adult biomass \n (kg km^-2)'
) +
ggtitle("Prediction (fixed effects + all random effects)",
subtitle = paste("maximum estimated biomass density (whole time series) =", round(max(exp(p$est))/1000))
)+facet_wrap(~Year,ncol=2)
## Warning: Removed 1169 rows containing missing values (geom_raster).
